Dynamical heterogeneities in a supercooled Lennard-Jones liquid 

Walter Kob 1 , Claudio Donati 2 , Steven J. Plimpton 3 , Peter H. Poole 4 , and Sharon C. 

Glotzer 2 

1 Institut fur Physik, Johannes Gutenberg- Universitdt, Staudinger Weg 7, D-55099 Mainz, 

Germany 

2 Center for Theoretical and Computational Materials Science, and Polymers Division, National 
Institute of Standards and Technology, Gaithersburg, Maryland, USA 20899 
3 Parallel Computational Sciences Department, Sandia National Laboratory, Albuquerque, NM 

87185-1111 

^Department of Applied Mathematics, University of Western Ontario, London, Ontario 

N6A 5B7, Canada 

(February 7, 2008) 
We present the results of a large scale molecular dynamics computer sim- 
ulation study in which we investigate whether a supercooled Lennard-Jones 
liquid exhibits dynamical heterogeneities. We evaluate the non-Gaussian pa- 
rameter for the self part of the van Hove correlation function and use it to 
identify "mobile" particles. We find that these particles form clusters whose 
size grows with decreasing temperature. We also find that the relaxation time 
of the mobile particles is significantly shorter than that of the bulk, and that 

this difference increases with decreasing temperature. 
PACS numbers: 61.43.Fs, 61.20.Lc , 02.70.Ns, 64.70.Pf 

Recent NMR experiments have shown that the relaxation in supercooled liquids is not 
homogeneous, i.e. that there are regions in space in which the relaxation of the particles is 
significantly faster (or slower) than the average relaxation of the system Subsequently, 
this result has been supported by optical spectroscopy, forced Rayleigh scattering and fur- 
ther NMR experiments [0. However, these types of experiments are unable to determine 
the nature of these "dynamical heterogeneities," and consequently details such as size are 
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unknown. Extrapolations of probe size sensitivity to heterogeneous dynamics indicate a 
typical heterogeneity size on the order of 2 - 5 nm in the vicinity of T g |3[] . 

Dynamical heterogeneities have also been observed in computer simulations [|j] . However, 
these simulations were restricted to two dimensions and since it is expected that the dynamics 
of particles in two and three dimensions is significantly different, it is not clear whether 
the dynamical heterogeneities observed in these simulations have a counterpart in three 
dimensions. By analyzing the trajectories of monomers in a Monte Carlo simulation of a 
dense (d = 3) polymer melt, Heuer and Okun || showed that in this system dynamical 
heterogeneities occur on short length scales, but the nature of the heterogeneities was not 
explored in detail. Thus despite the experimental evidence for the existence of dynamical 
heterogeneities, their microscopic properties are unknown and thus phenomenological models 
are often used to interpret experimental results |J. In this Letter, we study a simple, glass- 
forming liquid to investigate whether dynamical heterogeneities can be observed in a 3-d 
system and, if so, to determine their properties. 

We investigate a binary (80:20) mixture of 8000 Lennard- Jones particles consisting of two 
species of particles, A and B. The interaction between two particles of type a, (5 e {A, B} 
is given by V af3 (r) = 4e a/3 [(a a/3 /r) 12 - {(r a p/rf] with e AA = 1.0, a AA = 1.0, e AB = 1.5, 
cab = 0.8, tBB = 0.5, and <jbb = 0.88, with a cuttoff radius of 2.5a a p. Note that the 
AB interaction is stronger than both the AA and BB interactions, a fact which will be 
important in the subsequent discussion of the results. We report all quantities in reduced 
units, i.e. length in units of a AA , temperature T in units of e AA /kB, and time t in units 
of J & AA m I e AA , where m is the mass of either an A or B particle. We study the system 
at 10 different values of T ranging between 0.550 and 0.451. At each T, the system was 
equilibrated for a time longer than the a-relaxation time before evaluating the quantities 
presented below. At the lowest T, quantities were evaluated for over 4 x 10 6 time steps. All 
simulations were carried out in the microcanonical ensemble. More details on the simulation 
can be found in Ref. [0. 

The dynamics of this model has been characterized in detail in previous simulations 
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performed at different temperatures and at constant density ||. In particular, it was found 
that at low T the dynamics is described well by mode-coupling theory || with a critical 
temperature T c « 0.435 and a critical pressure P c ~ 3.03. In the present work we approach 
the point (T C ,P C ) via a different path than that used in Ref. [Q, on a straight line in the 
T — P plane along which density increases with decreasing T [|IU| . It has been shown |7|] that 



along this path of approach to the critical point, the behavior of the relaxation dynamics 
is very similar to that found along the constant-density path of the previous simulation |§ , 
providing evidence that the thermodynamic path of approach to the critical point does not 
significantly change the nature of the divergence of the relaxation time, and hence how the 
system vitrifies. Hence, we expect the results presented here to be independent of the details 
of the approach to the glass transition, and so in the following we use T alone to characterize 
the different state points. 

To detect the presence of dynamical heterogeneities, we investigate the time dependence 
of the self part G s (r,t) of the van Hove correlation function JTTJ] for the A particles, where 
r is the distance traveled by a particle in a time t. To a first approximation G s (r,t) has 
a Gaussian form but deviations from this form at intermediate times have been observed 
in simulations of glass forming liquids P,[T2|Jr3| and are thought to reflect the presence of 



dynamical heterogeneities ||14|| . Such deviations can be characterized by the non-Gaussian 
parameter of G s (r,t), a 2 {t) = 3(r 4 (t))/5(r 2 (t)) 2 — 1 |I5| . Fig. [I] shows the time dependence 



of a 2 for the A particles at three different temperatures. We find that: (i) on the time 
scale at which the motion of the particles is ballistic, a 2 is zero; (ii) upon entering the 
time scale of the /3-relaxation a 2 starts to increase; and (iii) on the time scale of the a- 
relaxation, a 2 decreases to its long time limit, zero. We observe that the maximum value of 
a 2 increases with decreasing T, which is evidence that the dynamics of the liquid becomes 
more heterogeneous with decreasing T. Furthermore, we find that the time t* at which this 
maximum is attained also increases with decreasing T. 

To determine the reason for the strong increase of a 2 in the /3-relaxation regime, we 
compare G s (r,t) with the distribution that is obtained from the Gaussian approximation, 



i.e. by assuming that G s (r,t) is given by G 9 (r,t) = (3/27r(r 2 (t))) 3 / 2 exp(— 3r 2 /2(r 2 (t))), 
where (r 2 (t)) is the mean squared displacement of the particles. In Fig. ^| we show (G s (r, t) — 
G 9 (r,t))/G 9 (r,t) for t = t*, where t* depends on T (see Fig. 0). For small and intermediate 
values of r (r < 0.6) the relative difference between G 9 S and G s is less than a factor of three. 
However, for larger r, G 9 S underestimates G s significantly. The discrepancy increases strongly 
with decreasing T in that the normalized difference becomes as large as 10 8 at the lowest 
T (see inset of Fig. Q). Thus we find that in the supercooled liquid there is a significant 
number of particles that have moved farther than would be expected from the Gaussian 
approximation [12j| . We define r* as the larger value of r such that G s (r*,t*) = G 9 (r*,t*), 
i.e. r* is the value of r at which the normalized difference starts to become positive and very 
large (see Fig. [|). We thereby define "mobile particles" as A particles that have moved 
farther than a distance r* within a time t* . With this definition, the total number of 
mobile particles at any T studied is a few hundred (out of 6400 A particles), and thus 
constitute approximately 5% of the system. We note that we find the results presented 
below concerning the properties of the mobile particles to be relatively insensitive to the 
details of this definition. 

Snapshots of the configuration of the mobile particles show that these particles tend 
to form clusters, i.e. they are not randomly distributed throughout the system. The spa- 
tial correlation between mobile particles is shown in Fig. |3], where we compare (cf. inset) 
9AmAm( r ) and gAA{r), the radial distribution functions for the mobile particles and for the 
bulk, respectively. (In the following, "bulk" refers to all of the A particles.) We find that at 
short and intermediate distances (r < 4) the mobile particles are more strongly correlated 
than the bulk. This is demonstrated more clearly by computing the ratio gAmAmij') / 9aa{j'), 
which is shown in Fig. |3| for three different T. From this figure we see that with decreasing 
T the relative correlation between the mobile particles increases; i.e., the relative height of 
the first nearest neighbor peak (r m 1) increases quickly and the ratio decays more slowly 
as a function of r if T is decreased. At the lowest T, the size of the cluster is on the order 
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of 3 - ^gaa [0- If we assume a molecule of diameter 0.4 - 0.5 nm we find that the clusters 
have a size of about 1 nm, which is in rough agreement with experimental expectations |||3| . 
We note that at small wave-vectors the partial structure factors for the bulk do not show 
any indication for the presence of these clusters. Thus it is perhaps not surprising that no 
evidence for the presence of such clusters was found from the structure factors measured in 
the neutron scattering experiments of Leheny, et al. ||17|| . 

What is the effect of these clusters of mobile particles on the bulk relaxation dynam- 
ics? To explore this question we compute the incoherent intermediate scattering function 
F( Am \q,t) for the mobile particles and compare it with that for the bulk particles, Fj A \ 
These correlation functions, shown in Fig. f| for three T, are calculated at a wave-vector 
q = 7.2, which coincides with the location of the main peak in the structure factor ||. 
From this figure we see that _Fj Am - ) decays faster than F^ A ' and that the ratio between the 
relaxation time of the two correlation functions increases with decreasing temperature. (The 
relaxation time could, e.g., be defined as the time it takes a correlation function to decay to 
e _1 of its initial value.) This ratio is approximately 3 for the highest T and approximately 
10 for the lowest T. It is not unreasonable to extrapolate that at temperatures close to 
T g the ratio between the relaxation times become as large as 10 2 — 10 4 , similar to values 
reported from experiments It is also interesting to note that the a- relaxation time 

of Fj" 4 ™) is on the order of the end of the /3-relaxation of the bulk. This suggests that the 
relaxation of these clusters might be related to the /3-relaxation of the bulk. 

Also shown in Fig. ^ is the probability P(t) that a particle which was mobile at time t = 
is still mobile at time t, at the lowest T investigated (bold dotted curve). We define P(t) 
by ((N(t) - N(0) 2 /N A )/(N(0) - N(0) 2 /N A ), where N A is the total number of A particles 
and N(t) is the number of particles that were mobile at time t — and still mobile at 
time t. A related correlation function for the least mobile particles has been measured in 
experiments We see that P decays on the time scale of the intermediate scattering 

function of the mobile particles, demonstrating that the lifetime of a cluster is of the order 
of the relaxation time of the particles which constitute the cluster. However, the lifetime is 



significantly shorter than the a-relaxation time of the bulk. More details on the dynamics 
of the particles within the clusters will be given elsewhere ||16|| . 

We find that the existence of clusters of mobile particles is related to small, local equi- 
librium fluctuations in composition |]IB|| . At all T, the pair correlation function gAmB^r) 
between the B particles and the mobile particles is smaller than the bulk quantity gAB{ r ) 
for r < 3. Thus mobile particles have fewer B particles in their vicinity than do generic A's. 
Because in this system the attractive interaction between A and B particles is stronger than 
either the AA or BB interaction, the presence of a B between two A's lowers the potential 
energy, giving rise to an effective attraction between the A's. A particles in a 5-rich region 
can thus be expected to have a reduced mobility. A particles in a B-poor region, however, 
will have a reduced effective attraction between them, resulting in a higher mobility [19| . 

Via this mechanism we expect that this sort of dynamical heterogeneity will occur in other 
fragile glass-forming systems, since local equilibrium fluctuations arising in the arrangements 
(packing) of the molecules will always be present . Energetically favorable packings will 
reduce the local mobility, while less favorable packings will enhance the local mobility. Such 
correlations have been observed in recent spin glass simulations ET] . 
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FIGURES 

FIG. 1. Non-gaussian parameter a.2 versus time t for T = 0.550, T = 0.480, and T = 0.451. The 

arrows mark the location of the maximum, i.e. of t*. 

FIG. 2. (G s (r,t)-G 9 s (r,t))/G 9 s (r,t) versus r for t = t* for T = 0.550, T = 0.480, and T = 0.451. 

The arrow marks the location of r* for T = 0.451. Inset: the same quantity on a logarithmic scale. 

FIG. 3. Inset: radial distribution function gAmAm{f) and gAA{f) for the mobile and bulk parti- 
cles, respectively, for T = 0.451. Main fig ure^ Ratio between (jAmAm 

(r) and gAA^f) for T = 0.550, 

T = 0.480, and T = 0.451. 

FIG. 4. The incoherent intermediate scattering function for the mobile (bold lines) and the 

bulk particles (thin lines) for T = 0.550 (dashed line), T = 0.480 (solid line), and T = 0.451 

(dashed-dotted line). Bold dotted line: Probability P that a particle which is mobile at time zero 

is also mobile at time t for T = 0.451. 
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